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Abstract 


The radiative heat transfer inside a low-density carbon fiber insulator is analyzed using a three-dimensional direct 
simulation model. A robust procedure is presented for the numerical calculation of the geometric configuration factor 
to compute the radiative energy exchange processes among the small discretized surface areas of the fibrous material. 
The methodology is applied to a polygonal mesh of a fibrous insulator obtained from three-dimensional microscale 
imaging of the real material. The anisotropic values of the radiative conductivity are calculated for that geometry. The 
results yield both directional and thermal dependence of the radiative conductivity. 
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Systems 


1. Introduction 


For highly porous media subjected to high temper- 
ature and high temperature gradients, radiative energy 
transport within the material is dominant. This radiative 
heat transfer can be expressed in the form of a radiative 
conductivity [1]. In such a model, the thermal conduc- 
tivity is highly dependent on the temperature. An ex- 
ample of an application where this type of heat transfer 
is important is in evaluating the design of Thermal Pro- 
tection Systems (TPS) for planetary probes and space 
vehicles [2, 3, 4]. 

In a previous study [5], a three-dimensional direct 
simulation model was developed to calculate the radia- 
tive conductivity of fibrous material using a simplified 
geometric model with properties resembling those of a 
carbon fiber insulator. It was shown that the radiative 
conductivity was a function of the temperature and the 
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orientation of the fibers because of the anisotropic na- 
ture of the test geometry. 

In the present paper, a three-dimensional model ob- 
tained from X-ray microtomography [6, 7] is used so 
that the geometric properties of the fibrous material at 
the microscopic scale can be better captured as manu- 
factured. More specifically, a rendering of FiberForm® 
is utilized. FiberForm® is a carbon fiber insulator that 
serves as the preform to certain TPS materials [8]. To 
demonstrate the validity of the approach documented in 
this paper, the effective conductivity tensor of the sam- 
ple is calculated, and compared to the experimental re- 
sults. 


2. Geometric configuration factor calculation and 
analysis 


The calculation of the geometric configuration factor 
(GCF) is an important step when determining the radia- 
tive exchange between two surface elements (or “sur- 
faces”). The computation of GCF involves integrating 
over the solid angles subtended by each surface ele- 
ment, if they are in sight to each other. In the absence 
of an exact analytical solution for the GCF, numerical 
integration is required. Since complex geometries re- 
quire a discretized geometry with many small surface 
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area of an internal enclosure [m7] 

area of a triangular surface element [m7] 
coordinate of the centroid [m] 
geometric configuration factor (GCF) 
geometric factor [m] 

length of the enclosure [m] 

number of surfaces within enclosure 
heat flux [W m7] 

distance between two surface elements [m] 
temperature [K] 

X,Y,Z spatial Cartesian coordinates [m] 


Greek Symbols 
radiative conductivity [W m7! K7"] 
principal-axis geometric factor [m] 
principal-axis radiative conductivity [W 
m7! Ky 
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elements to preserve accuracy, a robust procedure is re- 
quired. Theoretical details and general expressions of 
the GCF numerical calculation were presented and ver- 
ified in a previous study [5]. 

The numerical rendering of the surface of a complex 
geometry is often performed using a mesh of flat, tri- 
angular elements. In the limit of transparent gas, the 
energy exchange between the individual elements, con- 
sidered to be isothermal, opaque, diffuse and gray, con- 
tributes to the radiative energy exchange. 

To ensure that all radiative contributions are properly 
accounted for, the geometry is enclosed in a cubic fidu- 
cial volume (see Fig. 1(a)). Each boundary wall, called 
internal enclosure (IE), is composed of the same mate- 
rial as the fibrous material (isothermal, opaque, diffuse 
and gray). 


2.1. Surface-to-surface GCF calculation 


The GCF is calculated for three different cases: (1) 
triangle-to-triangle, (2) IE-to-triangle and (3) IE-to-IE. 
The expression corresponding to each case is listed in 
Table 1, where 6; is the angle between the normal of 
the surface elements and the line of length S;; joining 
the centroid of the triangles i and j, of area a; and aj, 
respectively. Each IE surface is divided into a set of 
uniform subsurfaces {AA} (see Fig. 1(b)) to perform the 
discretized integration. The size of AA is taken to be 


Stefan-Boltzmann constant [W m~? K~*] 
angle between the normal and S$ 
emissivity 
error due to grid resolution 

Superscripts 

eff effective 

sol solide 


Subscripts 

€,m,n directional indices 

III opposite side enclosure walls indices 
l, J,k surface index 

IP In-plane 

TTT through-the-thickness 

irr irradiation 

rds radiosity 

fin final 

ini initial 








(a) Cubic enclosure contain- 
ing the complex geometry, 
with internal enclosure (IE) 
surfaces identified by num- 
bers A. 


(b) Square mesh of subsur- 
face AA superposed on an 
internal enclosure (IE) sur- 
face of area A so that }AA = 


Figure 1: Schematic figure of the specific experimental geometry em- 
ployed in the model. 


small enough so that the results are considered mesh 
independent, with a relative error of less than 10~’. The 
reciprocal relations corresponding to each case is also 
listed in Table 1. 


2.2. Acceleration of the GCF calculation 


For N surfaces, owing to the reciprocity relation, only 
N(N — 1)/2 GCFs need to be calculated by integration. 
The calculation can be further accelerated by taking ad- 
vantage of the geometric properties of the surface: 


Table 1: Geometric configuration factor and reciprocity relations 
used in the model. Case (1) Triangle i to triangle 7. Case (2) Subsur- 
face k of IE / to triangle 7. Case (3) Subsurface k of IE 7 to subsurface 
€ of IE j. 


Case Geometric configuration factor Reciprocity 


cos 6; cos 6; 
1 F.3 =4;—{, — AP jj = ajF jai 
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AA cos 6; cos & 
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Field-of-view. A non-zero GCF exists if and only if sur- 
faces i and j are in sight to each other. Defining 6 as 
the angle between the normal of the surface element 
and the line joining the centroid of two triangles, then 
Fi; = F jj = Oif 0 ¢ [O, 2/2]. 


Obstruction. Incase of isotropic emission (diffuse-gray 
approximation), radiation exchange between two sur- 
faces is possible if and only if their line of sight is not 
obstructed by an opaque object. This occurs when an 
Opaque surface intersects the straight path from cen- 
troid 7 to centroid j. If this condition is met, then 
Fj a F jj = 0. 

It is worth noting that these two procedures do not 
necessarily guarantee that the radiation is completely 
blocked. In some cases based on the size and orientation 
of the triangles, the energy can be partially exchanged, 
since portions of surfaces are still in sight of each other. 
This inaccuracy may be mitigated by using a fine mesh 
of good quality [9], with triangles of approximately uni- 
form area. 


2.3. Refinement of the GCF calculation 


As the GCF is calculated, two additional procedures 
are performed to ensure that the discretization of the do- 
main is sufficiently refined: 


Bounding. If the calculated Fj_,; is greater than 1, then 
the mesh is considered to be too coarse. In the first equa- 
tion of Table 1, since cos 6; cos 6; is always € [0, I], it 
is possible to ensure that the mesh is sufficiently refined 
by requiring that a ;/(aS =) < | over all triangles. 


Closure. In an enclosure, the sum of all view fac- 
tors Fj_,; from a given surface i must equal to | (ie. 
Dj i+; = 1). Due to the numerical errors associated 


with the discretized surface elements, the summation is 
usually smaller than | (ie. >); Fi; = 1 — & where 
é; ~ 10°%). To remediate this, the value of &; is dis- 
tributed among all configuration factors of i such that 
the summation becomes 1. 


3. Radiation exchange and effective heat conduction 
model 


In a black-body radiation approximation, where ab- 
sorptivity and emissivity (€) coefficients are equal, the 
balance equation between radiosity g;a, and irradiation 
Gir, at temperature T for the i” surface element is given 
by [10]: 

Ards,i = coT; +(1- E)Girr,i ) (1) 


where o = 5.67 x 10°° W/(m’-K*) is the Stefan- 
Boltzmann constant. Quantity gir; 1s determined by 
summing over the portions of radiation within the fidu- 
cial volume that reaches surface 7. Using the reciprocity 
relation, Eq. (1) becomes: 


Gris = €0T} +(1—8) ) FinjQrisj. 2) 
J 


For the IE, a re-radiating boundary condition is applied 
such that g;qs = Girr, Which can be written as the follow- 
ing: 


rds,i = >, Fi jas, j . (3) 
J 


Equations (2) and (3) form a N x N system of equations, 
which can be solved to obtain the values of radiosities 
{rds,i} . 

In order to determine the effective conductivity of a 
volume due to radiative transfer, the energy flux pass- 
ing through the fiducial volume is modeled assuming a 
heat conductivity type relation. As illustrated in Fig. 3, 
radiative conductivity xk is expected to show a direc- 
tional dependency. Therefore, the law of heat conduc- 
tion (Fourier’s law) for an anisotropic medium is used: 


Gn = De KmnOnl, m,n € {x,y, Zz}. (4) 


n 


Here, 0, denotes derivative operator (i.e. 0/0,, 0/0, and 
O/0,). The Kmn elements are determined using a system- 
atic procedure. A linear temperature gradient ranging 
from T7; to T; is applied along the n—direction, such that 
the i” triangular surface is assigned temperature: 


Cui 


T; = T7; — AT L, (5) 





Here, AT = Ty; — T; > 0, C,,; and L, are n-component 
of i" triangular centroid coordinate and the distance be- 
tween the two enclosure walls in the n—direction, re- 
spectively. Also, a re-radiating boundary condition is 
imposed for the enclosure walls normal to the n (m # n). 
This temperature gradient is only imposed in the n— 
direction, therefore leading to Onz,T = O. The pla- 
nar averaged heat flux vector g,, for the imposed tem- 
perature is calculated using the difference between the 
radiosity from one enclosure wall, and the irradiation 
flux from the opposite wall. Therefore, from Eq. (4), 
for each n—direction, three conductivity elements Ky, 
(where m = x, y, and z) at T = (7, + T7;)/2 can be 
determined as: 


Kmn(T) = —Gn(T) 





Ly 
AT 
Based on the Rosseland formulation, the radiative 
conductivity is expected to exhibit a cubic polynomial 
dependence on temperature [10, Ch. 12]. Therefore, the 
Kmn elements are fitted to a cubic polynomial of temper- 
ature, 
Knn(L) = KnnéoT”, (6) 


where K is the geometric factor associated with the 
Rosseland model. 

It is important to point out that in this methodology, 
only the Rosseland temperature-dependency formula- 
tion is used, and not the approximate model [11]. 


4. Results 


4.1. Geometry 


The fiducial volume is obtained from a subset of a 
micro-CT scan of FiberForm®, obtained during a pre- 
vious study [7, 6]. The fibrous volume is segmented into 
a surface mesh of 56 268 triangular elements, illustrated 
in Fig. 2. properties of Fig. 2. 


4.2. Geometric configuration factor 


For illustrative purposes, a binned “3D map” of the 
bilateral GCFs of the geometry is shown in Fig. 3. Each 
value is obtained by averaging the 562687 F;_, j ele- 
ments over 1007 intervals. For instance, F;_, j Values for 
(i, 7) = (0,0) — (561,561) are averaged and mapped to 
the bin (0,0) and the Fj; values for (7, j) = (562,0) > 
(1 123,561) are averaged and mapped to bin (1, 0), and 
so on. The heterogeneous structure of FiberForm® 
manifests itself through the non-uniform distribution of 
the GCFs which consequently leads to a non-uniform 
energy distribution over all the fiducial volume. 





Figure 2: (Color online) Three-dimensional microscopic rendering of 
FiberForm® extracted from X-ray computed microtomography. The 
volume is contained in a 0.1 x 0.1 x 0.1 mm? cube, and the fibers have 
an average diameter of 10.6 um. 


4.3. Radiative conductivity 

The analysis for calculating the radiative conductiv- 
ity (Section 3) is performed for AT = 1.0 K over the 
model for a temperature range of Tini < T < Thin , where 
Tini = 300.0 K and 7%, = 4000.0 K. For the emissivity, 
a uniform value of ¢ = 0.85 is used [12]. The calcu- 
lated temperature-dependent conductivities are shown 
in Fig. 4. Because the conductivity tensor is known to 
be symmetric [13], the averaged temperature-dependent 
values of both k,,,, and K,,, are used to calculate the off- 
diagonal elements. As expected from the Rosseland 
model, these results follow the cubic temperature de- 
pendency of Eq. 6. Therefore, the geometric factors 
Kym associated with the best fit obtained from a cubic 
curve are extracted and presented in Table 2(a). The 
calculated values illustrate the anisotropic behavior of 
the selected volume since the diagonal elements are all 
different. 

The radiative conductivity tensor possesses three mu- 
tually orthogonal principal axes (i.e. the orientation of 
the geometry where the off-diagonal terms are zero) 
which are associated with three real eigenvalues. The 
Rosseland factor A for the principal-axis conductivity 
elements J are listed in Table 2(b) so that 


An(T) = AmeoT?. (7) 
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Figure 3: (Color online) Calculated geometric configuration factors 
have been plotted in a bin-averaged over 562 x 562 consecutive (i, j) 
combinations. 


Table 2: (a) Geometric factor reconstructed from Rosseland model, 
and (b) geometric factor reconstructed from Rosseland model along 
the principal axis. 


(a) K [10° m] 
xX y Z 
x 2.072 0.207 0.130 
y 0.207 1.526 0.436 
z 0.130 0.436 1.856 
(b) A [10-7 m] 
1 2 3 
2.347 1.213 1.894 


4.4. Discussion 


The values of the thermal conductivity calculated us- 
ing the present method can be compared to the values 
provided by the manufacturer [14], as well as those 
available in the literature [8]. In these two references, 
the values in the through-the-thickness (TTT) direc- 
tion is presented. The TTT direction is defined as 
Kmm < Knn> Kee, Where Knn © Kee 1S the in-plane (JP) di- 
rection. For FiberForm®, these directions are the results 
of the orientation of the fibers due to the compression 
axis. Since these results are obtained experimentally, 
they include all the physical phenomena that contribute 
to conductivity: the radiative conductivity, the gas phase 
conductivity (argon), and the solid conductivity. For the 
radiative conductivity factors calculated in Table 2, the 
TTT direction corresponds to the smallest value, the y— 
axis. 

To properly compare the two values, the radiative 
conductivity calculated in this work is added to the solid 
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Figure 4: (Color online) Radiative conductivity matrix elements as a 
function of temperature. 
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Figure 5: (Color online) Comparison of the total effective conductivity 
with experimental results presented by Tran et al. [8] and Fiber Mate- 
rials Inc. [14]. The total effective conductivity includes the radiative 
conductivity calculated in this work, and the solid conductivity om 


calculated by Panerai et al. [7] 


conductivity value in the TTT direction calculated by 


Panerai et al. [7], so that «KS... = krrr + K%,. To 
calculate «%%.,, they use a 1.0 mm?* cubic sample of 


FiberForm® in air, at 273.0 K. Even though it is ex- 
pected that the solid conductivity slightly decreases with 
an increase in temperature [15], a constant value of 
5 = 0.203 W/(m-K) is used here. In addition, cor- 
responding eigenvalue A = 0.165 W/(m-K) is cal- 
culated from the radiative conductivity tensor presented 
in [7], so that A =Arrrt ae is used for compari- 
son. 

The comparison is presented in Fig. 5. As can be 
seen, the manufacturer and experimental data, over the 
entire range of temperatures, are in very good agree- 
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Figure 6: (Color online) Ratio of the solid and radiative conductivity 
as a function of temperature: the red region illustrates the temperature 
range where both terms are within 10% of each other, indicating that 
they both should be included in the effective conductivity calculation. 


ment with the values bounded by the total effective con- 
ductivity Krrr + oe and corresponding principal axis 
total conductivity Arrr+ AS Experimental data in the 
other directions are not readily available, but a factor of 
Kip/Krrr = 2.4 1s mentioned in Ref. [8], and a factor of 
2.0 in Ref. [16]. Comparing to these ranges of values 
would also give a reasonable agreement for the two k7p, 
corresponding to the k,, and x,, values calculated here. 

Figure 6 presents the xrrr to es ratio as a function 
of temperature, highlighting the respective contribution 
of the two terms. It can be seen that for temperatures 
lower than 650 K, the radiative conductivity is likely 
negligible, accounting for only 10% of the effective con- 
ductivity. Likewise, at temperatures higher than 3000 K, 
radiative effects are dominant and the solid conductivity 
can be neglected. 

Finally, it is important to point out that the 
temperature-dependent gas phase conductivity and ra- 
diative absorption are not accounted for in the present 
work. Moreover, the geometry of the sample used for 
the radiative conductivity, Krrr, 1s different from the one 


used in Panerai et al. [7] to calculate a 


5. Conclusion 


The use of a X-ray microtomography generated 
geometry for computing radiative conductivity was 
demonstrated. A robust algorithm to calculate the ge- 
ometric configuration factors among the planar surface 
elements of test geometry using triangular mesh gen- 
erated by microtomography was constructed. The cal- 
culated values highlighted the heterogeneous behav- 
ior of fibrous material. The net radiation exchange 


model combined with Fourier’s law was used to inves- 
tigate the radiative energy flux distribution within the 
FiberForm® micro-structure. The developed model ac- 
counts for both the angular and thermal dependence of 
the radiative conductivity. The extracted values of the 
radiative conductivity tensor elements compare remark- 
ably well with the experimental data, and manifest the 
anisotropic behavior of the material. The results also 
show that the cubic polynomial formulation of Rosse- 
land model is valid for this case. 
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